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FOREWORD 


This report discusses the numerical solutions obtained for the simplest 
three-dimensional configuration, a flat plate, using the simplified Navier- 
Stokes equations. A separate report Is being written to present the results 
for an axial corner formed by two perpendicular plates. The leading-edge 
problem, which Is essentially two-dimensional, has been studied using both 
the exact ana the simplified Navler-Stokes equations, and will be described 
In another report. These program listings and user's guide are to be docu- 
mented in the fourth report. 
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COMPUTATIONS OF SUPERSONIC VISCOUS FLOW OVEF 
A FINITE-WIDTH PLATE 

ABSTRACT 


Finite-difference methods are applied to solve the parabolic Navier-Stokes 
equations for the flow over a finite-width plate at 0° and 10° angles of 
attack. The methods are developed on the basis of the operator-factorization 
concept resulting in the split of a three-dimensional equation Into successive 
two-dimensional equations. Two versions of the factori zation , backward and 
centered Implicit schemes, are used and their results are compared. Available 
numerical solutions and experimental data obtained at low Reynolds number con- 
ditions are also used for comparison. The backward implicit method provides 
a more successful solution, which ranges from the merged-layer to the strong- 
interaction regimes. The present study also reveals the detailed structures 
of the shear layer around and behind the side edge. 
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1 . INTRODUCTION 


The problem of calculating viscous flow over simple configurations has been 
of practical Importance in fluid dynamics In connection, for example, with 
the design of high-speed flying vehicles and with the investigation of flow 
phenomena in test facilities. It is known that the viscous and heat conduct- 
ing effects not only dominate the flow within a narrow layer adjacent to 
the solid boundary, but also play an important role in shaping the flow struc- 
tures near the leading edge, around the side edge, and behind the object. 

These flow fields, being far more complex and interesting than their two- 
dimensional counterparts, are not amenable to most of the existing methods of 
analysis due to some basic difficulties. In the first place, the governing 
equations, known as the Navier-Stokes equations, although providing an excel- 
lent description of the flow field, are an elliptic type for which the numeri- 
cal techniques are lengthy and cumbersome. Secondly, because of the occur- 
rence of shocks and shear layers resulting in steep gradients of flow proper- 
ties, sufficiently fine resolution is required to maintain acceptable accuracy 
as well as a stable solution. Thus, even in some situations in which numeri- 
cal solution is feasible, the computation cost is prohibitively high for fre- 
quent applications. Consequently, there have been some attempts to solve the 
flow problems with slightly different governing equations which possess para- 
bolicity in one of the independent variables, usually along the coordinate 
defining the main flow velocity. This set of equations, termed as the para- 
bolic Navier-Stokes equations, is to be integrated spatially and requires sim- 
pler and efficient numerical techniques, compared to that for the elliptic 
Navier-Stokes equations, the object of the present work is to develop a new 
algorithm to solve the parabolic equations and to evaluate its validity using 
a problem of the flow passing a finite-width plate. 

There exists strong evidence that the parabolic Navier-Stokes equations 
can be used successfully to describe a large class of flow problems in which 
the physical dissipations and the pressure gradient can be neglected in the 
direction of the main flow. For example, ref. 1 and 2 have presented solu- 
tions of the incompressible flow inside a rectangular duct, and ref. 3 and 4 
have discussed their techniques for solving the compressible flow over 
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a cone, an axial corner, and a finite-width plate. Although the problem 
studied may be different In flow property and In geometry, the governing equa- 
tions employed have a close resemblance In the underlying assumption that the 
streamwlse second-order derivatives can be omitted from the complete Navler- 
Stokes equations. Furthermore, *he numerical techniques developed are typi- 
fied by finite-difference integrations marching from one plane to the next, 
both normal to the streamwlse coordinate. The implicit schemes (ref. 1,2,4) 
are preferable over the explicit ones, owing to their unconditional stability 
which allows a larger increment In the marching direction regardless of the 
fact that steep gradients necessitate fine discretization of the space on the 
integration plane. The combination usage of iterative and implicit schemes 
(ref. 3) Is particularly appealing when the second-order derivatives are only 
present In one direction. In this technique and the ADI technique (ref. 4), 
the solution of a sparse-tanded linear system of equations resulting from 
the direct application of implicit schemes to the governing equation Is re- 
placed by the solution of a tridiagonal system. The advantages achieved in 
more efficient computation far outweighs the extra work involved. This also 
suggests a further exploitation of the method of splitting (ref 5) to con- 
struct other versions of difference techniques that may have more desirable 
features in the stability, accuracy and in the efficiency for practical compu- 
tations of three-dimensional flowfields. The essential concept of splitting 
is to seek the solution of a multidimensional equation in two or three steps, 
each involving the solution to a two-dimensional equation. In terms of dif- 
ference operators, this idea can be worked out readily by means of factoriza- 
tion, similar to that for an algebraic equation. 

The other part of the numerical algorithm involving the formulation of dif- 
ference equations and the solution to the nonlinear equations also receive 
considerable attention in this study. A mesh of cells, instead of points, is 
used because of its recognized compatibility to satisfy the conservative laws, 
both locally and globally (ref. 6). The cell formulation is recommended for 
flow containing large gradients and for separated flows. In order to solve 
the nonlinear difference equations, an iterative procedure is developed which 
retains the conservative property after the convergence is reached. A pre- 
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scribed tolerance on the order of the truncation error In the difference 
scheme is used to control the linearization error. These aspects are as im- 
portant as the difference schemes to the success of computation, and hence 
deserve equal attention in the development of a numerical algorithm. 

As a test case for the new algorithm, the flow past a finite-width plate at 
0° and 10° angles of attack is considered. This problem was solved in ref. 4 
using an ADI technique and its solution was among one of the earlier success- 
ful applications of the parabolic Navier-Stokes equations. However, it appears 
from the published results that the computation was terminated prematurely at 
a station quite close to the leading edge and that the inviscid-viscous interac 
tions near the side edge are not as strong as anticipated. These difficulties 
may be caused in part by the fact that the ADI scheme is only marginally stable 
and therefore not suited to deal with solutions containing rapid variations, 
and in part by the use of a uniform mesh system which would fail to detect any 
sharp change in properties around the side edge. 

The following discussion is divided into five sections, each concentrating on 
one aspect of the numerical solution. They are given In the order of the 
finite-difference methods, nonlinear difference equations, governing equations 
and boundary conditions, mesh system and stability, and the flowfield results. 
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2. FINITE DIFFERENCE METHODS 


Consider a model partial differential equation In the conservative-law form 

F x ♦ Q y ♦ H z « ° (2-D 

where F, G, and H are column vectors and also the function of a column vector 
V. x, y and z are the Cartesian coordinates. 


If the backward implicit scheme Is applied directly to eq. (2-1), the resulting 
difference equation Is 



.1 Al+l 

Jk ' Ay^ vJ+Sk 




( 2 - 2 ) 


Subscripts j and k denote the locations of the cell point, j+Sk and jk+S denote 
the location of the cell line. Superscripts 1 and i+1 refer to the previous 
and present values. The computation space defined in the y-z plane Is discre- 
tized to have a mesh of cells, each having variable Ay^ and Az^. The integra- 
tion increment, Ax i , is also determined prior to solving eq. (2-2). Due to 
the nonlinear form of eq. (2-2), a 1 i neari zation procedure Is adapted to ren- 
der a linear system of equation in terms of V^. This system of equations has 
a r oarse-banded coefficient matrix, and requires lengthy computation. In prac- 
tice, yr, ADI or an Iterative-Implicit technique in which the solution of a tri- 
diagonal coefficient matrix is employed. The other alternative technique is 
the method of fractional steps, whereby the solution to eq . ( 2-2 ) is obtained 
from two equations 


F* ■ F 1 
h jk h jk 

Ax.j 

( s w ■ 

■ ‘i-i.i 

* 0(4xJ) 


F^ » F* 
jk jk 

Ax.. 

("Ik ' 

■ >sd 

*°H) 

(2-3) 


The two successive steps constitute a cycle of the calculation, which can be 

efficiently performed just as the ADI technique since k and j are fixed in the 

first and the second equation, respectively. ,o 
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A simple analysis can be carried out to demonstrate the equivalence between 

o 

eqs.(2-2) and (2-3) to the order of (Ax^) . Assuming F-G-H-V, eq.(2-2) can 

be written In the operator notation such as 


LV 1 * 1 . V 1 


(2-4) 


and eq.(2-3) becomes 


L y V* » V 1 , l 2 V 1+1 ■ V* 


(2-5) 


where 


LV - -aVj. sk * «» J+W - BV^ ♦ BV Jk+li 


V ' -“Vw * V * ° V J*,k 


L k V ■ - 8 V Jk-S * V jk + 0 V Jk+S 


with a = Ax./ Ay^ and 6 * Ax^/Az k 


By substituting the second Into the first equation of eq.(2-5), It leads to 


L y L z V - L z L y V - LV ♦ 0 


K) 


Second-order accurate schemes may be constructed in the same manner, 
example, using the centered scheme, the difference equations are 

f Jk ■ F jk • W- ( G W • G W + G W - G j-Sk) * °K) 

Sr (v. - "jk-s + H ji^ • <->.) + °K) 


For 


pi +1 

jk 


F 3 


(2-6) 
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or, In the operator notation 


L y V* - 1/ , y,™ - L Z V* (2-7) 

where 


and 

r i v ■ • ? ^ 

Eq . ( 2-6 ) can be shown tc be equivalent to the difference equation of eq . ( 2-1 ) 
using the centered scheme In ora step. Eq . ( 2-6) is mo.'e accurate than eq.(2-3) 
for smooth values of V. However eq.(2-3) is preferable when V has discontinui- 
ties because the backward scheme provides numericil dissipation to smooth out 
unwanted small wave-length osci lotions. 

The validity of approximating a multidimensional difference operator by a 
sequence of two-dimensional operators has not received full theoretical 
justification it the operators a*e nonlinear. Nevertheless, this method has 
been used for solutions of various nonlinear problems and has also been used 
successfully in the numerical applications. Finally the nonlinear operators 
do not commute with each other In general, the following equations are used 
instead. Eqs . ( 2- 3 ) and (2-6) In the actual applications become 


l j( l/ L z v'* 2 . v< 

(2-8) 

L z( L y ) 2 ^ Vl * 2 ■ L z( L y ) 2 V’ 

(2-9) 


to maintain the desired formal order of accuracy. 
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3. SOl-riON TO THE NONLINEAR DIFFERENCE EQUATIONS 


Eq.(2-1) and Its difference equation, eq.(2-6), are nonlinear, so a method of 
linearization Is devised to reduce the equations Into a soluble form. It may 
be advantageous to linearize the difference equations directly as shown below 
because the conservative-law form can be retained throughout the Iterations. 

Assume there exist certain relationships between the column vectors F, G, H 
and V, V y , V^ t where the subscripts denote the partial derivatives. 

dF « PAdV 

dG • PBdV - PCdVy - PDdV^ 
dH - PEdV - PFdV z - PGdVy 

where P, A, B, C, D, E, F and G are square matrices, PA can be Interpreted as 
the Jacobian matrix of F with respect to V. etc. Making use of these relations, 
Iteration formulas are obtained according to the Newton-Raphson technique. 

V 4+1 > v 4 ♦ 6V 

F 4+1 - F 4 ♦ (pa) 4 6V 

G 4+1 « G* + (Pb) 4 6V - (PC )*■ 6V y - (PD) 4 6 V z 

H 4+1 « H 4 + (pe) 4 6V - (pf) 4 6V z - (PC) 4 6V y (3-1) 

Here, l indicates the iteration count. Updated functions are given on the left 
side of equations, while previous functions are on the right side. Since new 
Jacobian matrices are computed at each iteration, the convergence of eq.(3-l) 
is quadratic. In general, no more than five iterations are needed to ensure 
that 6V's be less than a reasonable magnitude of tolerance e. The prescribed 
value should be greater than the round-off error, but smaller than the trunca- 
tion error due to the finite-difference approximation of eq . ( 2- 1 ) . A check is 
made to see if | 6 V | <e after each iteration. If every point on a fixed j th or 
k th line satisfies the condition, this line will be dropped out from the itera- 
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tive process. The iterations are terminated when every line has converged 
solutions. This procedure is also known as the Guass-Seidel line iteration 
technique as the updated value is used immediately in the iteration procedure*. 

The correction vector 6V is determined by solving a linear system of algebraic 
equations, which is obtained after substituting eq.(3-l) into the first equa- 
tion of eq.(2-6). The system of equations has the following form 


b 2 6V 2 + c 2 V 


3 



a j 6 v j-i + y v j + c j^ v j+i 

a o 6v j-i + b j 6v j = d a 



(3-2) 


with j=3,4...J-l. The coefficients in the first and the last equations have 
incorporated the boundary conditions, hence they are slightly different from 
the coefficients in the second equation. A general expression of these coef- 
ficients can be derived as a function of the Jacobian matrices and the mesh 
spacings. 


a, = 


b o 


“iMjV 

- MS + (“2 ■ e 2 )( PB )j - ( pc W( y i 

‘ ( PC )j+!/( y j+l ' y ’ 


' y j-l) 


j) 


Ay 4_ If'': - F 1 ') + G 1 .., - G 1 . ,+ G 1 ,, - G* 
X i VJ 3/ 3+h 3-H 3+k J- 


d j " 2"Ax ■ 


i-% 


(3-3) 


Note that these groups of coefficients eq.(3-3) and eq.(3-2) are used to solve 
for <SVj at a fixed k th line. The mixed-derivatives involving values across 
the kth line do not enter into the expressions. A similar type of equations 
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and coefficients car. be derived for 6V k at a fixed J th line. These equations 
constitute a block trldlagonal system for which an efficient algorithm (described 
In ref.101 Is readily adopted for solution. It Is of Interest that as 1 6V ^ | -►t , 
|dj must approach to e at the same rate. The recovery of the first equation 
eq.(2-6) at the end of Iterations Implies that the conservative laws are fol- 
lowed strictly. The weighting parameters are a^«l/(l+R), cyc^R, R-AY^/AYj ; 
likewise and $2 are obtained using R»AYj + j/AYj. 

Substituting eq.(3-I) Into the ‘;’1rst equation ot eq . ( 2-3 ) results In a block, 
tridiagonal system of linear equations identical to that shown In eq . ( 3-2) . 
However, the coeffirlents differ somewnat from those given In eq . ( 3-3 ) . bj 
and dj have the following re atlons Instead. 

b j ■ • sf K *K • •»)<«)} - MjV (»i - »j.i) 

- ( pc )jv( y j*l ■ y j) 

d j • § ft ■ F i) + • g U < 3 -*> 

Thereft-e the secc..d-order accuracy in Ax can be gained in eq . ( 3-2 ) and eq . (3-3) 
with little extra work. 

An alternate approach to linearize eq.(2-9) is to solve for V directlv from 

0 

V . The Newton-Raphson technique leads to relations such as 

f^ 1 *• f* ♦ (pa)*(v £+1 - V £ ) 

and the resultant coefficients closely resemble those in eq.(3-3) and (3-‘0, 
but tne accuracy may not be as good because the round-off error is larger duo 
to V | >> ( 6V | . Furthermore, the derivation of coefficients at the boundaries 
becomes more complex (ref. 6). (5V 1 * <5Vj=0 will be replaced by ' v j s Vj • 
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4. GOVERNING tQUATIONS AND BOUNDARY CONDITIONS 


The Integration of the governing equations Is performed along the x-axIs in 
the y-z plane. The computational space is defined by 0<y<H and 0<z<B , in 
which the nlate lies on the z-axis and 0<z<W. 


The functions in eq . ( 2-1 ) 

F x + G y + H z " 0 (2--) 

are defined by 


"pu 

G = 

"v 

pu^+p 


PVU+Tyx 

puv 


PV 2 +o y 

| puw 


P VW+T yz 

|(pe+p)u_ 


[i )e+o y ) u+T yx u+T y z + q y _ 


H * 

"pw “1 


I pwv+x zy 
pw 2+o z 

(pe+o z )w+T zx u+T zy v+q z 


(4-1) 


The boundary conditions are given as 

y=0, u*Au y , v=w a 0, e*e w - ^ jj- e y 

0<z<w 

w<z<h u y a w y = e y = p y = p y * 0 , v=0 (a=0) , v y =0(a/0) 

z=0, u z = v 2 ■ e z = p z = p z = 0 , w=0 

y=H or u*U cos a, v*-U sin a, w=0, e=e , p=p , p*p (4-2) 

* oo oo oo' ooao 
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where p, u, v, w, p and e are, respecti vely , the density, velocity components 

In (x,y,z), pressure and th® Internal energy. The total Internal energy Is 

2 2 2 

e»e+ 0.5(u +v +w ). p indicates the pressure to be treated differently from 
those appearing in the stress components, which are defined as follows. 

yx y’ yz zy i y z) zx z 

o y ■ p - (X ♦ 2y) v y - Aw z 

o 2 « p - (A + 2y) w z - Av y (4-3) 

The heat fluxes are obtained from the expressions: 



where 


2A + 3y * 0, k = C p p/P r , e * C y T and y * C p /C v 


u is the viscosity coefficient given for air by 


V 


2.27xl0” 8 


t 3/2 

T+19876 


lb/ sec 
■ — 

ft*' 


(4-5) 


k is the coefficient of thermal conductivity and is related to the Prandtl 

number', P = 0.71 for air. C and C are the specific heats at constant pres- 
r p v 

sure and volume. Their ratio, y, is equal to 1.4 for perfect air. T denotes 
the temperature, a is the angle of attack. 


Eq . ( 2-1 ) is suppletiented by the equation of state 

p = (y-1 )pe (4-6) 
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The boundary conditions nave taken Into account the slip In velocity parallel 
to the wall and the jump in temperature. A, the mean free path of air, can 
be estimated by 


A 




U 


(4-7) 


Two other Important parameters used very often are the Reynolds number. Re, 
and the viscous-inviscid Interaction parameter , They are defined as 



where is tne f. 'ei-stream Mach number and C is the Chapman-Rubesin constant, 

C » uT Jv J . 

To facilitate the numerical computations, three new vectors are Introduced 
as 


p 

• V w * 

"p ” 

, v, = 

rp,n 


y 

y 

z 

z 

u 




u. 



y 


z 

V 


v„ 


V 



y 


z 

w 


w 


w 



y 


z 

e 


e , 


e. 

i— — J 


Ly_j 


L 2 J 


(4-8) 


The Jacobian matrices of F, G, and H with respect to V, V , and V are 

y z 


PB 


3F 



3V 



3G 

3V 

PC 

3G 

3H 

PF 

3H 

3V 

* 3V. 


•sSSSS* 


OY 


- 13 - 


0 V 

ge/p 0 
0 0 

ave 0 


ygpe 0 



0 0 0 

que P+9P 0 


c ■ p 0 
0 


F .I 0 

0 


0 ygu/p 


ge/p 0 
gwe 0 


ygpe gpw 


X+2y 0 
0 ygy/P r 


2 2 2.2 a . Y _i = 0, and * 0. 

and q * u + + w . g - Y- * . a p be 

. in the derivation of linear system of equation 

These expressions are use The Jacob1a „ s assoc1a{e d with mixed 

*" d also in the / re not needed, so they are omitted from 

derivatives, i.e., gy- and aV" > a 

. z y 

eq. (4-9) . 



5. MESH SYSTEM AND STABILITY ANALYSIS 


The computational domain outlined In section 4 Is divided into a mesh of rec- 
tangular cells, each of which has variable length and height. The smallest 
cell is located at the leading edge, whereas the largest cell Is at the oppo- 
site corner of the region. Integers are used to designate the cell number; 
half-integers designate the line between two cells. This system of cells Is 
determined prior to the computation by means of a procedure to be described 
In the following. First, an estimate is made of Ay^ or Ay 2 using the smallest 
value of the boundary-layer thickness at X*L and the mean free path of the 
free stream, l.e., Ay^ ■ m1n(6,A), where 6 * 0. 31 / and A given In the 
previous section. Then an exponential function is employed to determine Ay^ , 
wnere j«3,4...Jp . This ^unction relates the transformed coordinate y in 
the computational space t") the physical coordinate, by j * (e c - 7 -l) h/(e c -l ) , 
0<y<l , 0<y£h. C is the parameter to be determined iteratively in order to 
satisfy the requirement that y^ *0.5 Ay 2 , y„ +1 * yj+O.SAy^ arid y Jp * h. Jp 
is the number of fine cells in the y direction. On top of this system of fine 
cells, there is a system of coarse cells. The height of those cells are uni- 
form and obtained from A y^ * (H-h)/ ( JL-JF) with y^ +1 * y^+O.SAy and y JL * H, 
where JL is the total number of cells along the j-axis. 


The width of cells is nonuniform in the region near the side edge and uniform 
elsewhere. The procedure to determine az r and z k are similar to that for Ayj 
and y^ except that for this case Az k$ » Az k$ _ 1 are specified to insure that 
fine resolution is provided at the side edge, ks designates the selected 
location of the side edge. Likewise, A*, and x^ are estimated using the re- 
quirement that Ax-| * A^ , and Xj *IL where ILis the total number of integration 
steps. 


The mesh of cells is constructed to properly resolve the gradients of flow pro- 
perties that may occur near the leading edge and around the side edge. It is 
unnecessary to be concerned with the efficiency and the stability of the numeri- 
cal computations because as for an implicit scheme, these aspects are no longer 
dictated by the CFL condition controlling Ax^. A simple analysis of stability 
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is presented to show that the present algorithm indeed is unconditionally 
stable. Consider a set of linear equations corresponding to eq.(2-l) 


AV+BV-CV ■ 0 

* y yy 


(5-1) 


A, B, and C are given in section 4. 


Eq . ( 5-1 ) is composed of the hyperbolic und the parabolic equations, and their 
difference forms obtained from the centered scheme are 


' cr • -I) 

*; w ■ -;) 


Ax 


.i+4 _ V 1+S\ 


4yJ B j V j-W 


(5-2) 


Ax 


1 


<v 


c ] ( v j*i - v r s * v !-t) <5 - 3) 


The amplification matrices for these equations are obtained using the Von 
Neumann's method (ref. 7) as 


r _ I - \/*ia A ^ B s1n 2 y 
G h • . i jr - 

I + v *ia A B sin y 

G , I - /mb a - 1 C sin 2 J 
I + /mb A’ 1 C si n 2 J 


(5-4) 

(5-5) 


where 



Ax. 




Y = JAyj and I is the unit matrix. 


The eigenvectors are obtained by solving polynomials of the following equations. 


U H 

-A-’ 

B| 

* 0 

I* 

-A* 1 

C 

* 0 

p 


1 



origin M. 
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Thus 



(5-6) 

(5-7) 


since tney are real variables, that implies |G^|<1 and |G p |<l. Thus eqs.(5-2) 
and (5-3) are always stable. 


If the backward scheme Is used, the numerator In eqs.(5-4) and (5-5) are re- 
placed by an identity matrix. Following the same procedure, the same conclu- 
sion is arrived at. 


It is worth mentioning that numerical instability could arise if the stream 
wise pressure is treated on the same basis as other variables in eq.(2-l). 
For this case, F and A become 




2 2 2 

They will have complex components if u +c <a , wherea is the local sonic 
speed. Further, and in eq.(5-9) are also complex. Therefore eqs.(5-2) 
and (5-3) become unstable if p is used instead of p. This stability problem 
was known in earlier work of solving the parabolic N-S equations (ref.* ) and 
P x was eliminated. from the governing equations. A better approach is to employ 
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the pressure gradient obtained at the upstream station, but the Improvement 
of accuracy Is not noticeable If Ax^s are made nonuniform (ref. 8). 

The efficiency of Implicit schemes is not necessarily higher than that of ex- 
plicit schemes. The advantage of selecting the larger step Increments and 
extremely fine cell sizes Is penalized by a relatively large amount of compu- 
tation time solving the nonlinear equations. In fact, more Iterations are re- 
quired In the procedure described In section 3 for regions where properties 
change rapidly, compared to regions near the undisturbed free stream. In 
general, the efficiency of a particular scheme depends upon the nature of the 
problem to be solved, and varies appreciably with the equation formulation and 
the performance of the nonlinear solver. For a leading-edge problem, ref. 4 
indicates that the ADI scheme Is substantially more economic to use than the 
straightforward explicit scheme. The present algorithm should be comparable 
with the ADI in the aspect of efficiency; however, there is no evidence avail- 
able to verify it at this time. 
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6. FLOW F I ELD OVER A FINITE-WIDTH PLATE 


The numerical computations have focused on a case that was studied and reported 
In ref. 4. A limited amount of data was also available from a test facility 
consisting of a nozzle placed Inside a Mach 12 blowdown tunnel. Due to the 
large drop In stagnation pressure, extremely low densities are achieved in 
the test section. The free stream conditions upstream of the plate have a 
Reynolds number. Re * 300/in and a Mach number, M^ ■ 5.15. The free stream 
temperature Is T w ■ 230°R, while the wall temperature is T^ ■ 460°R. The 
numerical results were obtained from an ADI scheme using a mesh of points with 
constant spaclngs in the transverse and lateral directions , also In the stream- 
wise directions. Comparisons with the experimental data were made at X2l0; 
a station displays strong merged-layer characteristics . 

The computational domain and the mesh of cells were devised In a different 
manner than those used in ref. 4 In an attempt to properly represent the 
flowfield structures around the side edge. The size of the domain is given 
as follows: L«2 In., H-2.4 In., h*0.48 In., w*2.16 in. and b*3.03 In. The 
cell system is designated by the notation (JF)JLxKL(KS) , as (12)32x32(20), 
where JL and KL are the total number of cells along y and z coordinates. The 
numbers in the parentheses indicate the total number of nonuniform cells along 
y, and the location of the side edge. The resulting cell dimensions are 
Ay 1 »Ay 2 *0.024in. , and Az 20 =Az ig *0.0123 in. A schematic of the cell system is 
shown in fig. 2. The streamwlse increments are also nonuniform and determined 
by a procedure described in the previous section after assigning Ax^*0.024 in. 
and K=30, the number of steps. Note the smallest cell used is greater than the 
mean free path, which is A^O.03 in. for this case. 


Comparisons are first made among the theoretical predictions of the pressure 
distributions across the stream in the lateral direction. Figure 3 shows 
that there are significant differences in pressure values near the side edge, 
but a good agreement toward the plane of symmetry and the opposite side in the 
free stream. The centered solution indicates a very sharp drop of pressure 
right at the side edge; the backward solution predicts a moderate variation, 
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whereas the previous ADI results display a slight change in pressure across 
the side edge. It is very difficult to assess the accuracy yielded by these 
computations without data from other sources. By observation from the pub- 
lished results in ref. 4, the shear layer characteristics might have been par- 
tially ignored because the mesh system is relatively coarse, even though the 
computational domain is less than half of that used in the present computa- 
tions. On the other hand, the appreciable differences between the centered 
and the backward solutions are affected to some extent by the numerical dissi- 
pations of the schemes. The centered scheme is known to provide more accurate 
results if the solution is smooth, but may exaggerate the results for a rapidly 
varying solution. Further, due to its lack of damping, small -wave length dis- 
turbances tend to grow and eventually become outbound. Indeed, it was found 
in the present calculation that the centered solution encountered some diffi- 
culties with the iterative procedure after reaching x =6 *6 at K=25 when negative 
densities appeared in some cells neighboring the side edge. In contrast, the 
results obtained from the dissipative backward scheme are stable throughout 
the integration and appear to yield reasonable accuracy. For this reason, the 
backward scheme has been employed In the entire study. 

Other than the differences existing between theoretical pressure results, there 
is also slight disagreements with the available experimental data in the flow 
flux at the midspan of the plate. The reason for the deviation within the 
boundary layer is not yet determined. The theoretical predictions of flux agree 
with each other very well, except near the region where the maximum value occurs. 
These curves are not asymptotic to zero at y=0, since the slip velocity is used 
there. 

Figure 4 presents a series of pressure distributions vs. z-coordinate. The 
pressure starts out with the free stream value, slowly builds up its magnitude 
as the slip velocity reduces, and finally reaches the peak value p=5.25 p^ at 
X=19.5. At the beginning there is a smooth transition from the surface pressure 
to the free stream pressure outside of the plate. It is seen In figure 4(a) 
that the pressure distribution resembles closely that within a normal shock. 
Figures 4(b) through 4(d) show the detailed variations of pressure as the flow 
moves downstream. A discontinuity is developed within the shear layer sur- 
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rounding the edge, then dissipates and diffuses sideways as the shear layer 
grows thicker and reduces the magnitude of gradients. At the end of the com- 
putation, X-5.35 or X«2 In., the side edge has started to lose Its Influence 
on the local flowfleld. The mechanism producing the pressure discontinuities 
must originate from the strong interactions between the flow and the- plate. 

Due to the presence of a plate In the flow, the kinetic energy transforms 
directly Into the Internal energy on top of the plate. This process doe*; not 
take place Immediately because the rarefaction effect predominates at the 
leading edge. The flow over the plate adjusts Itself rather quickly, but 
slowly In the shear layer around the edge. Therefore, the pressure retains 
a higher value for a longer distance outside the plate until the physical 
dissipation fl.ially smooths out the nozzles. Unfortunately, no experimental data 
are available In this region, so confirmation of the prediction cannot be made 
at this time. Note that the side edge effect has almost reached the center- 
line as indicated In figure 4(e). This points out the possibility that the flow 
cennot be considered two-dimensional at X»L when the aspect ratio, W/L, Is less 
thar; unity. 

2 

The distributions of the coefficients of skin friction, t yM /SpU , and of 

w 

heat transfer, (q+uT^J/pJJ^fHj. -H^), are giver, in figures 5(a) and 5(b), re- 
spectively. As explained earlier, the shear layer Increases Its extent with 
decreasing values of x» and moreover the shear layer creates relatively higher 
heating at the edge compared to that at the center. Judging from these figures, 
however, the flow field remains two-dimensional at the centerline. 

Figure 6 shows the velocity vector, constant contours of w, pressure, and tem- 
perature plotted on the plane normal to the x-coordinate at x*5.35. It indi- 
cates an upward motion of flow above the plate. Near the edge, however, a side- 
ways motion predominates. Away from the side region, there are no transverse 
and lateral veloc.ty components. The lateral velocity is plotted In constant 
contours to show that the highest value occurs somewhere out in the free stream 
and that w changes rapidly at the side edge. The pressure contours exhibit a 
leading-edge oblique shock above and a side-edge shock next to the side edge. 

The temperature contours disclose the location of a distontinuity near the 
edge, which has a stronger intensity than the leading-edge shock. 
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More computations were carried out at angles of attack, a»10° and >10°, which 
represent the wind and the lee sides of the plate, respectively. The primary 
purpose Is to check the capability of the numerical algorithms under more 
severe conditions. It would be of Interest from the practical point of view 
to consider the entire flowfleld on both sides of the plate, but this neces- 
sitates extensive modification of the program while It adds little value to 
our objective. In the present approach, these two cases are treated Indepen- 
dently, with different boundary conditions Imposed at y-0 and b<z<w to approxi- 
mate the flow vector continuation. Other than this modification, free stream 
velocity becomes u a> * u ae COSci and v^«-Us1na. The computational region and the 
cell system are the same as those for a»0°. Despite the fact that this approach 
sacrifices physical velocity, the predictive potential of the algorithm can 
still be evaluated and assessed. 

Figure 7 shows the pressure distributions vs. z-coordlnate given at selected 
x-statlons. The wiggles In pressure distributions are believed to be caused 
by Insufficient space resolution, since flow properties change more drastically 
around the side edge. Regardless of the local numerical problems, the compu- 
tation was completed at x«2" for both cases. The x-component of shear stress 
and the heat fluxes behave In the same manner as those predicted for a»0°, 
hence they are not repeated. One additional feature, which bears no counter- 
part In the a»0° case. Is the flow separation across the lee side of the plate 
and a stronger lateral flow on the wind side. Shown In fig. 8 are the skin 
friction coefficients, , vs. z on both sides of the plate. The In- 

fluence of the side edge is felt at an earlier station on the center plane. 

It Is indicative that for a plate of aspect ratio equal to unity, the flow on 
the plane of symmetry is definitely not two-dimensional when the plate is 
placed at an angle to the free stream. The significance of the results ob- 
tained Is in the dScovery of separated flow immediately over the side edge. 

The lateral extent of separation is seen to be proportional to the strearrwise 
coordinate. The separation becomes stronger as x increases. Figure 9 shows 
that the transverse extent of the cross flow separation also Increases with x. 
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Figure 10(a) shows the lee side cross flow velocity vector obtained at x-2" and 
4". It exhibits a stronger lateral reversed flow region farther away from the 
leading edge. But, a complete vortex has not yet been established, probably 
due to the small angle of attack, which Is not able to produce a full vortex. 
Figure 10(b) Is the constant contour plots of lee side pressure at two stations. 
The pressure discontinuity adjacent to the side grows as does the lateral ve- 
locity. 

All the computations were performed on an IBM 370-168 computer using approxi- 
mately four minutes of CPU. The computation time required can vary from one 
case to another, because even with the same number of cells and streamwlse 
steps, some may need more Iterations in solving the nonlinear difference equa- 
tions than others. An understanding of flow characteristics is helpful in the 
selection of computational region and mesh, which could reduce the requirement 
of computation time and increase the accuracy of the results. 



7. CONCLUSIONS 


It has been shown that the finite-difference methods constructed using the 
method of fractional steps has promising potential In solving the parabolic 
Navler-Stokes equations. This study has concentrated on the computation of 
viscous flow passing a rlnlte-wldth plate, and utilized buth centered and 
backward difference schemes to Integrate the equations along the streamwise 
coordinate. The applicability of the centered «*.h»me Is found to be suscep- 
tible to shocks and shear layers In the flow field and hence not used in 
calculations. The solutions obtained from the backward scheme appear to be 
stable and reasonably accurate. The numerical algorithm is featured with 
the conservative-law formulation In the difference equations and with a lin- 
earization procedure that is also consistent with the conservative law. Re- 
marks pertaining to the flat-plate problem at flow Incidence are summarized 
as follows: 

1. This work has revealed the existence of a shear layer around the side 
edge of the plate. Its Intensity reduces from the peak value near the 
leading edge as the flow moves downstream, but It remains dlscernable 
well into the weak-interaction regime. 

2. The pressure distribution along the center line is more sensitive than 
other measurable quantities to the flow expansion about the edge. For 
the conditions used, the flow should be treated three-dimensional ly when 
the aspect ratio, W/L, is close to unity. 

3. A cross flow separation is predicted on the lee side at c s -10 J ; however, 
the complete formation of a vortex has not been established. 

4. The solution to this problem can be obtained from the elliptical Navier- 
Stokes eauations using a time-dependent method, but it would require 100 
times as much amputation time. 

5. The present method of solution restricted to a class of viscous flow 
problems which do not involve a local reversed flow region in the direction 
of the main stream and the streanwise pressure gradient is negligible. 
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